Freezing of 4 He and its liquid-solid interface from Density Functional Theory 
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We show that, at high densities, fully variational solutions of solid-like type can be obtained 
from a density functional formalism originally designed for liquid 4 He. Motivated by this finding, 
we propose an extension of the method that accurately describes the solid phase and the freezing 
transition of liquid 4 He at zero temperature. The density profile of the interface between liquid 
and the (0001) surface of the 4 He crystal is also investigated, and its surface energy evaluated. The 
interfacial tension is found to be in semiquantitative agreement with experiments and with other 
microscopic calculations. This opens the possibility to use unbiased DF methods to study highly 
non- homogeneous systems, like 4 He interacting with strongly attractive impurities/substrates, or 
the nucleation of the solid phase in the metastable liquid. 

PACS numbers: 68.08.De, 64.70.Dv, 64.30.+t, 67.80.-s 
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I. INTRODUCTION 

Helium crystals represent a fascinating system in which 
general properties of crystalline surfaces such as the equi- 
librium crystal shape, surface phase transitions (rough- 
ening) and elementary mechanisms of the crystal growth 
can be studied over a wide temperature (T) range, in 
principle down to absolute zero. Helium crystals can be 
grown chemically and isotopically pure, and with very 
few lattice defects. Helium crystal growth is fast in com- 
parison to that of other solids and it is characterized 
by the unique -quantum- phenomenon of crystallization 
waves, i.e., melting- freezing waves which can easily prop- 
agate on the liquid-solid interface at low temperatures. 
This allows an accurate measurement of the surface stiff- 
ness 7 and surface tension a and of their anisotropy . 
The research on the surface of helium crystals has been 
recently reviewed Q. 

A microscopic approach to the study of solid 4 He 
and liquid-solid coexistence at low temperature needs a 
fully quantum approach. Accurate Path Integral Monte 
Carlo Q (PIMC) in the 5-30 K temperature range and 
Green Function Monte Carlo Q (GFMC) calculations on 
solid 4 He have been reported in the past. More recently, 
the liquid-solid interface has been adressed within a 
Variational Monte Carlo (VMC) approach using shadow 
wavef unctions 

Density Functional (DF) methods @ represent a useful 
computational tool to study the properties of quantum 
inhomogeneous fluids, especially for large systems where 
DF provides a good compromise between accuracy and 
computational cost. Indeed, the T = properties of 
inhomogeneous liquid 4 He can be accurately described 
within DF theory (DFT) by using the phcnomenologi- 
cal functional proposed in Ref. [7j and later improved in 
Ref. §, and similar ones. They have been widely used 
in a variety of problems involving inhomogeneous 4 He 



systems like, e.g., liquid- vapor interface 0,01, pure and 
doped clusters [3, 0] , layering and prewetting transitions 
in films alkali atom adsorption on the surface of liq- 
uid 4 He and droplets 0], vortices in 4 He clusters |ll|. 
etc. 



In view of its conceptual interest on the one hand, and 
of the potential applications on the other hand, we have 
undertaken a fully variational study of the liquid-solid 
transition and coexistence of 4 He at zero temperature 
within DFT. Previous attempts on 4 He [Tl HI 111 ll^ 
or quantum hard spheres [TtI Hsj have been reported. 
They adapted the techniques initially designed to ad- 
dress the liquid-solid transition in classical systems (for 
a review see, e.g., Ref. 13]). They use trial functions 
(sums of gaussians) to parametrize the solid density, 
and resort to a second-order expansion of the energy 
around the liq uid density (Ramakrishnan-Yussouff (RY) 
approach) [T3, HJ, Ha, Ha] , or to a modified weighted den- 
sity approximation exact to second order for small per- 
turbations of the liquid 0, 0]. These methods need 
the density-density static linear response function x 01 
the liquid as input. In helium, these approaches have 
been found to be, at most, a reasonable starti ng p oint to 
qualitatively describe the freezing transition 12]; more- 
over, some of these attempts always give the liquid Q 
or the solid 0,0] as the stable phase. Additional prob- 
lems arise in treating the liquid-solid interface, because 
it is found to be too narrow to allow the use of coarse 
graining methods Besides, the use of trial function 
density profiles makes all these approaches of limited ap- 
plicability for systems characterized by the presence of 
both solid and liquid phases. 
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II. THE DENSITY FUNCTIONAL FOR SOLID 

4 HE 



We attempt here a fully variational DF description of 
both liquid and solid phases on the same footing. Our 
starting point is a simplified version of the Orsay-Trento 
(OT) functional @, where the energy of the 4 He system 
is written as 



e\p]= ^l dr ^V~pY 



drdrV(r)p(r')VLj(|r- 



+ f/ dvp(v)p 2 (v) + | J drp(r)p 3 (r) 

The first term is the usual quantum kinetic energy term, 
the second term contains a Lennard-Jones He- He pair po- 
tential VljM screened at distances shorter than a char- 
acteristic length h. In the third and fourth terms, the 
weighted density p is the average of the density over a 
sphere of radius h. These terms account phenomenolog- 
ically for short range correlations. The parameters h, C2 
and C3 (in the original formulation, h — h) are fixed to re- 
produce the experimental density, energy per atom and 
compressibility for the liquid at zero temperature and 
pressure @. The original non-local kinetic energy term 
proportional to Vp(r) • V'p(r') @ has been dropped, be- 
cause, as one might expect, it leads to unavoidable in- 
stabilities in the solid phase. Its neglecting causes the 
static response function x of the liquid to be only quali- 
tatively described. We should mention the existence of an 
alternative parametrization of this term that still repro- 
duces the experimental static function, while being free 
from instabilities As the resulting DF stems from 

a somewhat different strategy than current DF's do, we 
have preferred to keep using an OT-like functional, in- 
troducing the changes mentioned below. 

The equilibrium density p(r) is obtained by minimiz- 
ing E[p] with respect to density variations, subject to 
the constraint of a constant number of He atoms N. 
This is achieved by evolving in the imaginary time a 
non-linear Schrodinger equation for the order parame- 
ter \&(r) = ■J p(r), where the Hamiltonian is given by 
H = -h 2 V 2 /(2m) + U[p\. The effective potential U is 
defined in terms of the variational derivative of the last 
three terms in the energy functional Eq. $T§ with respect 
to * 0. 

The calculations are performed in a periodically re- 
peated supercell containing a fixed number of 4 He atoms, 
and the starting configuration is a superposition of gaus- 
sian profiles with arbitrary but small width, placed in the 
positions of the hep structure, which is the experimental 
solid structure of 4 He. Interestingly, for average densities 
corresponding to the stable liquid phase, the equilibrium 
density obtained is homogeneous, whereas for higher den- 
sities -those roughly corresponding to the solid- a solid 
4 He structure, characterized by a periodic density distri- 
bution with helium atoms at the lattice sites of an hep 
crystal, is found to be the stable phase. 'Solid' means 



r'|) 
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FIG. 1: Surfaces of constant 4 He density (p 
a 4 He hep crystal, viewed along the c-axis. 



0.029 A" 3 ) of 



here a highly inhomogeneous configuration characterized 
by regions, called 'atoms' for brevity in the following, 
where the density is very large. Such atoms are only 
slightly overlapping with the density tails of neighboring 
atoms. As an example, we show in Fig. ^the calculated 
density profile for one such solid structure, computed at 
the experimental freezing density pj — 0.0287 A -3 . 

We have studied other ordered structures as well. An 
fee lattice is also found to be a stable phase at typical 
solid densities, its energy being only slightly higher than 
that of the hep phase (by 0.02 K per atom). A simple 
cubic structure, at the same densities, is found instead 
to be unstable towards the homogeneous (liquid) phase. 
We want to point out that the correct solid structure is 
found even if the initial gaussian superposition has not 
the correct parameters. For instance, we found that by 
starting with gaussians, placed on a hep lattice with a 
lattice constant twice the equilibrium one, but such that 
the total normalization gives the correct 4 He density, the 
equilibrium hep structure is eventually recovered at the 
end of the minimization process. 

In spite of this encouraging result, a close analysis of 
the density dependence of the total energy of the hep 
solid 4 He, calculated using Eq. (JIJ, shows that it largely 
deviates from the experimental result. The reason is an 
unphysical, large density pile-up in the core region of the 
4 He atoms, such that the resulting solid is too dense. 
To fix this unrealistic behavior, we need a mechanism 
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TABLE I: Freezing transition parameters (a = 2.556 A). 



Method 



pia 3 p s a 3 Pf (bar) 



l/p (A 3 ) 



FIG. 2: Solid lines: calculated solid and liquid EOS; dashed 
line: double tangent construction. Solid circles are the exper- 
imental data for the solid EOS [2l|: the experimental EOS 
for the liquid is not shown since it coincides, by construc- 
tion, with that of the OT functional, whose agreement with 
experiment is excellent 0. 



that makes energetically costly such density pile-up. The 
simplest one is to add a 'penalty' energy term to the 
functional, which has the following form: 



E p [p] = C / p(v)f\p(v)]dv 



(2) 



Here /[p(r)] is a 'switch' function which becomes appre- 
ciably different from zero only when the density is larger 
than a predefined value p m . One such function is 



/[p(r)] = l + tanh{/3[p(r)- Pm ]} 



(3) 



where C, (3 and p m are DF parameters. The effect of 
E p [p] is to add to the effective potential U[p] the term 



V p (r) = C{f[p(r)]+Pp(r)[l~Unh 2 [/3(p(r)-p m )}} . (4) 

If C > and when p(r) > p m , this term acts as a repul- 
sive barrier, which forbids extra pile-up of the density. 
We can thus use C, /3, and p m as adjustable parameters in 
such a way to get agreement with the experimental equa- 
tion of state (EOS) for the solid phase. Note that there 
is a large freedom in the choice of the penalty term since, 
by construction, it has no effect whatsoever on the liquid 
structure, even in a highly structured liquid, because in 
practice p m is a factor > 10 than typical liquid densities. 
We have found that the choice C = 0.1 Hartree, (3 = 40 
A 3 , and p m = 0.37 A" 3 [U, yields an EOS for the solid 
in good agreement with experiments plj . as shown in 
Fig. 2. 



III. 



Experiment [22] 0.434 0.479 25.3 

GFMC [4] 0.438 0.491 27.0 

VMC [5] 0.449 0.456 27.0 

RY [12] 0.459 0.515 41.5 

RY [16] 0.435 0.513 25.7 

Present work 0.437 0.490 26.1 



THE FREEZING TRANSITION AND THE 
LIQUID-SOLID INTERFACE 



Once the EOS of the homogeneous liquid and the EOS 
of the solid are determined, one can easily locate the 
freezing transition by means of a double-tangent Maxwell 
construction. We find a freezing pressure Pf = 26.1 bar 
and a chemical potential at coexistence /if = 0.69 K. Our 
calculated values for Pf and the solid and liquid densities 
at coexistence are reported in Table[|] together with those 
obtained by other authors. 

We now turn to the liquid-solid interface of 4 He, char- 
acterized by the interfacial tension a and width £. 4 He 
is an example of a system for which a has been di- 
rectly measured, thanks to the possibility of propagat- 
ing crystallization waves at its surface. The value of a is 
0.17 x 10~ 3 N/m 1], but f is not known experimentally. 

The He liquid-solid interface is an example of a sys- 
tem whose direct simulations by means of exact Monte 
Carlo methods can be computationally very expensive. 
Recently, it has been tackled by means of shadow wave- 
functions VMC El . and also by the rescaled DF we have 
referred to before [lg. We have used the modified DF, 
Eqs. (1 1 12ft . to study the liquid-solid interface. To de- 
scribe the coexisting phases, we start from a slab of hep 
crystal 4 He in contact with a region of liquid 4 He, and 
minimize the energy functional to find the equilibrium 
configuration, imposing that during the minimization the 
chemical potential is constant and equal to the value 
/if = 0.69 K found from the double-tangent construction 
illustrated in Fig. |2 We consider explicitly the interface 
between liquid 4 He and the (0001) surface of the crystal 
(basal plane). 

Fig. |3 shows a view of the equilibrium density config- 
uration at the liquid-solid interface. A related quantity, 
the 4 He density averaged over a plane parallel to the 
basal plane, is shown in Fig. 0J The thickness £ of the 
liquid-solid interface has been estimated by evaluating 
the 10-90% width of the curve interpolating the maxima 
in the density profile shown in Fig. 21 We find 9.3 A, 
in close agreement with the VMC results. Our findings 
are compared in Table II with other experimental and 
theoretical results. 

From the calculated density profiles we estimate the 
liquid-solid surface tension as a = (E— [if N+PfV)/(2A). 
Here V is the volume of the supercell, and A is the surface 
area of the exposed (0001) face. A factor 1/2 appears in 
the equation to account for the two free surfaces delimit- 
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FIG. 3: Liquid-solid interface shown by means of equal den- 
sity contour lines (drawn between p = 0.02 A -3 and p = 0.05 
A -3 ) in a plane perpendicular to the interface plane. Con- 
stant density surfaces (at p = 0.08 A -3 ) are also shown to 
identify the atoms in the solid slab. 



ing the solid film in our slab geometry. Our calculations, 
which are three dimensional in nature, have been done 
by using a supercell accomodating ~ 6 layers of solid 
4 He and a thick liquid layer in contact with it. We find 
that the convergence with respect to the amount of liq- 
uid present in the supercell is rather slow (to obtain a 
accurately, one should have two very wide liquid and solid 
regions in contact). Our estimated value is a — 0.1 x 10 -3 
N/m |23|, which is a rather good result, of quality similar 
tothcVMC value @. 

TABLE II: Liquid-solid interface parameters. 

Method a (10~ 3 N/m) g (A) 

Experiment [1] 0.17 
VMC [5j 0.25± 0.1 10-12.5 
DFT [16] 0.47 5.6 
Present work 0.1 9.3 



IV. CONCLUSIONS 




10 20 3Q 40 



FIG. 4: Planar-averaged 4 He density plotted along the direc- 
tion normal to the interface. The solid is on the right, the 
liquid on the left. 



We have shown that DF's currently used for liquid 4 He 
admit fully variational solutions of highly inhomogeneous 
type that may represent localized helium 'atoms'. With 
a simple modification, one such DF is found to accu- 
rately describe both liquid and solid phases of He and 
their coexistence. Our scheme does not rely on any spe- 
cific trial function density profile for the solid phase, but 
gives an unbiased description of both phases on the same 
footing. This broadens the range of applicability of cur- 
rent DF, permitting to study in an unbiased way highly 
non-homogeneous 4 He systems, like droplets doped with 
strongly attractive impurities, as for instance alkali ions 
|24|. or 4 He on strongly attractive substrates, such as 
graphite or carbon nanotubes [2{|. Another potentially 
useful application is the study of the nucleation of the 
solid phase in the metastable superfluid prll ]. 
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through the program HPC-Europa Transnational Access. 
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